
function [p q] = implicit_theta(dt,beta)
global Ux Uy  h  g  p  q  dx dy L Nx Ny nt;

dt2= 1/dt;

L = Constructmatrix( Nx, Ny, Ux, Uy, dx, dy ,h ,g );

I = eye(2*Nx*Ny);
% Periodic Boundary Conditions
sol(1:Nx*Ny,1) =p;
sol(Nx*Ny+1:2*Nx*Ny,1) =q;

for n=2:nt+1;   
    rhs = (dt2*I - (1-beta)*L)*sol;
    sol=(dt2*I+beta*L)\rhs;    
    p = sol(1:Nx*Ny,1);
    sol(Nx*Ny+1:2*Nx*Ny,1) ;
    
%     if rem(n,100)==0
%         refreshdata(h1,'caller') % Evaluate p in the function workspace
%         drawnow
%     end
    
end
display('Completed Successfully');